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ABSTRACT 

We study the linear and nonlinear evolution of the tearing instability on thin current sheets by means 
of two-dimensional numerical simulations, within the framework of compressible, resistive magneto¬ 
hydrodynamics. In particular we analyze the behavior of current sheets whose inverse aspect ratio 
scales with the Lundquist number S as S' -1 / 3 . This scaling has been recently recognized to yield the 
threshold separating fast, ideal reconnection, with an evolution and growth which are independent of 
S provided this is high enough, as it should be natural having the ideal case as a limit for S — 1 - oo. 
Our simulations confirm that the tearing instability growth rate can be as fast as 7 ~ 0.6 r^ -1 , where 
ta is the ideal Alfvenic time set by the macroscopic scales, for our least diffusive case with S = 10'. 
The expected instability dispersion relation and eigenmodes are also retrieved in the linear regime, 
for the values of S explored here. Moreover, in the nonlinear stage of the simulations we observe 
secondary events obeying the same critical scaling with S , here calculated on the local, much smaller 
lengths, leading to increasingly faster reconnection. These findings strongly support the idea that in 
a fully dynamic regime, as soon as current sheets develop, thin and reach this critical threshold in 
their aspect ratio, the tearing mode is able to trigger plasmoid formation and reconnection on the 
local (ideal) Alfvenic timescales, as required to explain the explosive flaring activity often observed in 
solar and astrophysical plasmas. 

Subject headings: plasmas - MHD - methods: numerical. 


1 . INTRODUCTION 

Magnetic reconnection is thought to be the primary 
mechanism providing fast energy release, readily chan¬ 
neled into heat and particle acceleration, in astrophys¬ 
ical and laboratory magnetically dominated plasmas. 
Within the macroscopic regime of resistive magneto¬ 
hydrodynamics (MHD), however, classical reconnection 
models predict timescales, in highly conducting plas¬ 
mas, which are too slow to explain bursty phenomena 
such as solar flares in the corona or tokamak disrup¬ 
tions. In particular, the Sweet-Parker model (hereafter 
SP) o f two-dimensional, steady, incompressible reconnec¬ 
tion (Sweet 1958| Parker 1957) predicts a reconnection 
rate M = v/ca ~ where v is the speed of the 

flow entering the reconnecting site, ca the Alfven veloc¬ 
ity based on the field far from the sheet and S = Lca/v is 
the Lundquist number for a given magnetic diffusivity 77 
(L is the current sheet length or breadth, identified with 
the macroscopic scale), which can be as high as S ~ 10 12 
in the solar corona, if simply due to collisional resistiv¬ 
ity. Such a rate is way too slow to explain any of the 
impulsive phenomena described above. 

As first demonstrated by means of 2D MHD simula¬ 
tions by 


ing, SP- 


Biskamp (1986) however, stationary reconnect- 
ike sites become unstable once the Lundquist 


number exceeds a critical value of order S ~ 10 4 , and 
are subject to fast tearing modes and plasmoid forma¬ 
tion when their aspect ratio L/a becomes large enough, 
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also increasing the local reconnection rate. Recent de¬ 
tailed linear analyses and simulations have confirmed 


taney et al. 20U9 Hhattacharjee et al. 

2009; Cassak 

et al. 2009 Huang & Bhattacharjee 2010 

Uzdensky et al. 


In particular, the SP current sheet, of inverse 

aspect ratio a/L ~ S' -1 / 2 , in the presence of the typ¬ 
ical inflow/outflow pattern characterizing steady recon¬ 
nection, was shown to be tearing unstable with growth 
rates "/ta ~ S 4 / 4 , where ta = L/ca- For a recent review 
on the latest theoretical works on 2D reconnection and 
secondary island (plasmoid) instabilitie s, from MHD to 


Hall regimes, see Cassak & Shay (2012). 

The existence of instabilities with growth rates scaling 
as a positive power of S poses severe conceptual prob¬ 
lems, since the ideal limit, corresponding to S -A 00 
would lead to infinitely fast instabilities, while it is well 
known that in ideal MHD recon nectio n is impossib le. 


This issue was resolved by Pucci & Velli (2014) (PV 
hereafter), who studied the stability of current sheets 
with generic inverse aspect ratios a/L ~ S~ a . The au¬ 
thors found showed that a critical exponent separates 
current sheets subject to slow instabilities, with growth 
rates scaling as a negative power of S, from the un¬ 
physical fast instabilities scaling as a positive power of 
S. Indeed, for a = 1/3, they found the growth rate of 
the fastest reconnecting mode to become independent of 
Lundquist number. They therefore conjectured that cur¬ 
rent sheets should not collapse to aspect ratios greater 
than this critical value, at which point the instability, 
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which they called the “ideal” tearing mode, leads to 
Lundquist-independent reconnection. For this aspect ra¬ 
tio, current sheets have a thickness up to 100 times larger 
than a typical SP reconnecting layer, and the instability 
developed X-points and plasmoids, thus preventing any 
collapse to the standard SP current sheet or any other 
steady configuration with a > 1/3. This novel “ideal” 
tearing instability is very fast, with an asymptotic growth 
rate jta — 0.62, and leads to the sudden formation of 
several plasmoids. In particular, PV found the relation 
kL ~ S' 1 / 6 , with k the fastest growing wave-vector along 
the current sheet. 

In the present work we investigate numerically, by 
means of (compressible) resistive 2D-MHD simulations, 
the linear and nonlinear stages of the tearing mode for 
a current sheet at the critical thickness a/L = S ~ 1 / 3 . 
Several initial configurations are tested, from the Harris 
sheet with fluid pressure balance, to the purely force- 
free case, with magnetic field rotation inside the current 
sheet, and also different values for the asymptotic plasma 
beta. 

The goal of this paper is, on the one hand, to retrieve 
all the known linear results and scalings, namely the 
expected instability dispersion relation and eigenmode 
structure (within the range of Lundquist numbers acces¬ 
sible to our simulations, that is up to S = 10 7 to limit 
the computation time), on the other hand to explore the 
nonlinear regime of the “ideal” tearing instability for the 
first time. Our simulations provide further proof of the 
existence of such an instability, which is expected to set 
in during current sheet collap se arising in any turbulent 
scenario of plasma dynamic s (Lourci ro et al. 2009| |Ser- 
vidio et aL|2009 Rappazzo fe Farker||2U13 ). 


The paper is structured as follows. In section [2] we 
describe the set of equations, the initial conditions, and 
our numerical setup. Section [3] is devoted to the numer¬ 
ical validation of the linear theory of PV. In section [4] 
we show the nonlinear results. Section [5] contains the 
discussions and conclusions. 


2 . NUMERICAL SETUP 


We integrate the compressible, resistive MHD equa¬ 
tions numerically, with the adiabatic index 7 = 5/3, in 
the form 


dp 

dt 


+ V • (pv) = 0 , 


( 1 ) 


^ + (v-V)v= Vvp+(VxB)xB], (2) 

dt p 


dT 

~dt 


+ ( v 'V)T = ( 7 —1) 


-(V • v)T + 


1 |V x B| 2 ’ 
S p 


(3) 


1 

— = V x (v x B) + —V 2 B, (4) 


where S is the Lundquist number defined above and 
other quantities retain their obvious meaning. Physical 
quantities are normalized using Alfvenic units, namely a 
characteristic length scale L , a characteristic density po, 
and a characteristic magnetic field strength Bq/^/A/k (the 
background values measured far from the current sheet). 
Velocities are then expressed in terms of the Alfven speed 
ca = Bq/ y/Anpo, time in terms of ta = L/ca, the fluid 


pressure in terms of Bq/Ait. Note that we are using 
the energy equation © written for the normalized tem¬ 
perature T = p/p, where we use as a reference value 
To = {m/k^)c 2 A . With the given normalizations the 
Lundquist number is basically the inverse of the mag¬ 
netic diffusivity, namely S = caL/p —> p~ l . 

The initial conditions at t = 0 for our two-dimensional 
simulations of the tearing instability are different forms 
of the Harris current sheet configuration in which the 
equilibrium magnetic Held varies only in the x direction, 
reaching an asymptotic magnitude B = 1 far from x = 0, 
the plasma density is uniform p = 1 , and the pressure 
(and temperature) p = T = P/2 far from the current 
sheet, localized around x = 0 , the asymptotic plasma 
beta being a given parameter. Two types of equilibrium 
are considered: in the first case, the classical Harris sheet, 
the field has only one component, the magnetic pressure 
gradient is balanced by a temperature enhancement in 
the current sheet itself (requiring a local plasma beta 
of order one, regardless of the value of the parameter /3), 
whereas in the second case the magnetic field is in a force- 
free equilibrium, p and T are taken constant everywhere, 
the condition B 2 = const being preserved by the fact 
that the magnetic field rotates across the sheet, so that 
there is a non-vanishing component B z ^ 0 inside the 
current sheet itself. Introducing a new parameter (, the 
(normalized) maximum amplitude of the z component 
of the magnetic field, both equilibria can be described 
writing: 


B = tanh 


o 


and for the fluid pressure is 



1-C 

2 


£ sech ( — ) z, 

(5) 

\aJ 

2 , _ . 


-sed^J, 

( 6 ) 


with £ = 0 for the first case with B z = 0 in pressure equi¬ 
librium (PE hereafter), and £ = 1 for the force-free equi¬ 
librium (FFE hereafter). Intermediate cases of mixed 
fluid/magnetic pressure equilibrium with 0 < C, < 1 can 
are also equilibria. 

With these normalizations, it is the thickness a of the 
current sheet that defines the growth rate of the tear¬ 
ing instability: in the incompressible linear analysis by 
PV it has been shown that, when a ~ S' -1 / 3 , for suf¬ 
ficiently high values of the Lundquist number (larger 
than 10 7 — 10 8 ) the growth rate of the instability 7 , 
measured in terms of the macroscopic Alfvenic time ta, 
becomes independent of the magnetic diffusivity and of 
the order unity. For the set of simulations shown be¬ 
low the current sheet thickness has been taken precisely 
a = regardless of the equilibrium model chosen, 

i.e. the adopted values of /3 and 

The compressible, resistive MHD equations (l|4) are 
solved in a rectangular domain [— L X ,L X ] x [0, L/]with 
resolution N x and N y respectively. In the ^-direction, 
in order to resolve the steep gradients inside the current 
sheet using a reasonable number of grid points, we limit 
our domain to a few times the current sheet thickness, 
i.e. we set L x = 20a: this is a good compromise between 
the high resolution required inside the current sheet and 
the need to have boundaries sufficiently far from the re¬ 
connecting region. Along the 7 -direction the length is 
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chosen in order to resolve for the fastest growing modes 
of the instability (see below). 

The tearing instability is characterized by the expo¬ 
nential growth of modes with wavelength larger than the 
current sheet thickness, that is with ka < 1 (fc being the 
mode wave-vector along y). The length L y along the cur¬ 
rent sheet is thus adapted to cover the range of unstable 
modes, namely we choose L y = mX, where m is the num¬ 
ber of wavelengths A = 27r/fc that we wish to simulate in 
our numerical box. Both lengths are chosen to decrease 
with S, 

L x = 20 a = 20S'~ 1 / 3 , L y = m\ = m — S^ 1 ^ 3 , (7) 

KCL 

where the value of ka ~ S' -1 / 6 is of order 0.1 for S = 10', 
from the linear analysis by PV. The linear analysis in 
the next section is performed by taking sheet lengths 
for which only one mode (m = 1), the most unstable 
one, is excited, whereas for the nonlinear simulations we 
will choose L y so that several unstable modes are in¬ 
dependently excited (typically m = 4), so to allow the 
subsequent mode-coupling and inverse cascade (i.e. the 
merging of plasmoids). 

In order to trigger the tearing instability, the equilib¬ 
rium configuration is modified at t = 0 with velocity per¬ 
turbations of amplitude e ~ 10~ 3 (the rather large value 
speeds up the evolution) and wave-vector k = k y, where 
k is the same quantity appearing in equation (|Tb, namely 
the wave-vector of the fastest growing mode selected for 
the analysis, as expected from the linear theory. Along 
the ^-direction these velocities are concentrated at the 
current sheet location and vanish far from the current 
sheet. Moreover, the v x component is taken to be odd 
across the reconnection layer, whereas the y-component 
v y , is obtained by imposing the perturbation velocity 
field to be incompressible. The analytical expressions 
for the perturbations are 


v x =etanh£e -i cos {ky+<pk), (8) 

v y = e[2t; tanh£—sech 2 £ )e _? sin(fcy+<pfc),(9) 

where ipk is a random phase (for each mode k) and £ = 
x S 1 / 2 . 

The numerical simulations are performed by integrat¬ 
ing equations (l]|4| with an MHD code developed by our 
group. Along the current sheet, where periodicity is as¬ 
sumed, spatial integration is performed by using pseudo- 
spectral methods, while in the 27-direction integration is 
performed by the use of a f ourth-order scheme based 
on compact finite-differences (Lele 1992). The bound¬ 
ary conditions in the non-periodic direction are treated 
with the method of projected characteristics jPoinsot fc 
Lele|ll992l IRoe fc Balsara||1996| |Del Zanna et al.||2UUl| 
Landi et al.||2005|), here assuming non-refiecting bound- 
ary conditions. Time integration is performed using a 
third-order Run ge-Kutta method. Details of the code 
are described in Landi et al. (2005). 

The resolution is adapted to the Lundquist number we 
use: for S = 10 5 and S = 10 6 we choose N x = 1024 and 
N y = 128, while for S = 10 7 the number of cells in the x 
direction is increased up to N x = 2048. In the periodic 
direction we use N y = 128 for the single-mode runs (in 
order to reduce the computational costs as many Simula¬ 



ka 



ka 


Fig. 1.— The instability dispersion relation (growth rate as a 
function of fc, normalized against a -1 = S 1 / 3 ) for different values of 
the asymptotic beta (top panel: ft = 0.1; bottom panel: /3 = 1.6), 
and Lundquist numbers (red color: S = 10 5 ; green color S = 10®; 
blue color: S = 10'). Solid lines are the theoretical expectations, 
symbols are for numerical results (crosses: PE; squares: FFE). 


tions are required to reproduce the instability dispersion 
relation curves), while we take N y = 256 in the nonlinear 
reference simulation. We have verified that this relatively 
low resolution along the periodic direction is adequate, 
due the extreme accuracy of Fourier methods and the 
rather smooth gradients observed in the y direction. In 
spite of the relatively high values of S, in addition to 
instability, the equilibrium diffuses on time-scales which 
although long compared to the instability, are still suffi¬ 
cient to affect linear evolution leading to slig htly u nder¬ 
esti mate the growth rates of linear modes (Landi et al. 
2008). To avoid this, in the single mode linear analysis 
described below, the diffusion term of the initial equilib¬ 
rium is subtracted on the rhs of the induction equation 
at all times. 


3. SINGLE MODE SIMULATIONS: LINEAR 
ANALYSIS 

A first set of simulations of the tearing instability in 
current sheets with a = S ~ 1 - /3 is performed to confirm 
th e expected linear behavior, i.e. the scalings reported 
Pucci & Velli (2014[), and in particular the instability 


dispersion relation as a function of the model parameters, 
here reported in figure IT] (solid lines). 

Even in the presence of the general equilibrium in equa¬ 
tions (5][6) with ( = B 0z /B 0 ^ 0, it is easy to show that 
the linear analysis of the instability is unchanged with 
respect to PV (in the incompressible limit p = 1 and 
assuming perturbations in the x — y plane alone with 
d z =0). Th e govern ing equations for the tearing mode 
are still ( |Furth et al.||1963| ) 

7 K - fc2 w) = ik[B 0y (b x - k 2 b x ) - B 0y b x ], (10) 
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Fig. 2. — Profiles across the current sheet of the tearing instability eigenmodes b x ( b y is automatically determined by the solenoidal 
constraint), v y , and v x . Top panel: analytical linear theory; bottom panel: numerical results. All profiles are normalized to their maximum 
value. 



7 b x = ikB 0y v x + S 1 {b x "-k 2 b x ), (11) 

where we have assumed that perturbations are factorized 
as oc f(x) exp(7 1 + ky) and the prime denotes derivation 
with respect to x. Thus, even when ( / 0 there is no 
coupling of modes with the B 0z component, and the PV 
results for initial equilibria with a = 5 1-1 / 3 should re¬ 
main unchanged. Moreover, no dependency on (3 is ex¬ 
pected, as neither this parameter enters the instability 
equations above, thus the theoretical dispersion relation 
curves only depend by the choice of S. 

In spite of the expectations commented above for an 
incompressible situation, the results of a numerical sim¬ 
ulation can deviate from the analytical case, due to the 
compressible regime, to differences in the treatment of 
boundary conditions (see below), and in general to dis¬ 
cretization errors and other numerical approximations. 
Therefore, we choose to test the two limits of our initial 
equilibria for the current sheet, namely PE (£ = 0) and 
FFE (£ = 1). Moreover, we investigate both cases with 
/3 < 1 (f3 = 0.1) and /? > 1 (f3 = 1.6b where /3 is the 
asymptotic plasma beta in equation ([6]). Finally, three 
different values of the Lundquist number are tested here, 
namely S = 10 5 , 10 6 , and 10 7 , for a total of 12 sets of 
simulations of the linear phase of the tearing instabil¬ 
ity, with the aim of reproducing the expected dispersion 
relations numerically , as shown in figure [I] 

As anticipated in the previous section, lor each value 
of ka we vary L y while always selecting a single mode 
m = 1. The growth rate of the instability is computed 
by measuring the x-averaged amplitude of the compo¬ 
nent b x of the perturbed magnetic held ( b x = 0 at the 
initial time). In order to better compare with theoretical 
expectations, the PV eigenmode analysis have been re¬ 
done here for a limited region across the current sheet of 
—20 < x/a < 20, precisely as in our simulations. How¬ 
ever, additional discrepancies are still expected, since the 
v x eigenmode is forced to vanish at the boundaries in the 
PV calculations, whereas in simulations we impose non- 
reflecting boundary conditions. 

The first thing to notice by inspecting the computed 
dispersion relations is that, as predicted by the classical 
linear theory, for each value of S the curves have a max¬ 


imum at a given k, the peak location decreasing in k as 
S increases. The growth rate of the instability (normal¬ 
ized to the inverse of the large-scale Alfven time ta) has 
peaks ranging from 7 « 0.5 for S = 10 5 to 7 > 0.6 for 
S = 10 7 . 

In general we find that the simulations with (3 = 1.6 
(bottom panel) are more precise in matching the analyt¬ 
ical results than those with (3 = 0.1 (top panel), since a 
large beta is a condition closer to incompressibility (for¬ 
mally corresponding to an infinite value for the sound 
speed). Moreover, we find that simulations of the FFE 
scenario (squares) yield higher and usually more accurate 
values of the growth rates as compared to those employ¬ 
ing the PE settings (crosses): this is probably due to the 
fact that the purely force-free equilibrium leads to in¬ 
trinsically less compressible fluctuations. Finally, rather 
large discrepancies are observed for small scales (large 
values of ka), especially in the PE case. 

In figure [2] we plot the profiles of the perturbations b x 
(by is determined by VB = 0), v y and v x , all normalized 
to their respective maximum, across the current sheet in 
the x direction. In the top panels we show the analytical 
results, that is the eigenmodes of the linear analysis (here 
the PV calculations have been recomputed by imposing 
v x = 0 for x = ±20a), and in the lower panels we re¬ 
port the numerical solutions for a simulation in the FFE 
scenario with S = 10 7 and ka = 0.10, at a given time 
of the linear evolution of the instability. In order to re¬ 
cover the theoretical eigenmodes, velocity and magnetic 
field perturbations are shown with a 7r/2 shift in ky, as 
expected. Notice the steep gradients arising within the 
current sheet (|x| < a), where a high resolution is needed 
to resolve the small scales developed during the instabil¬ 
ity evolution. 

As seen, the eigenmodes are very well reproduced: the 
magnetic field perturbations are identical to the analyti¬ 
cal expected ones, while in the velocity perturbations the 
only major difference is, as anticipated, due to the non- 
reflecting free-outflow boundary conditions, that do not 
force v x = 0 at x = ±a and result in a slightly different 
profile even in the vicinity of the reconnecting region. 
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Fig. 3.— Nonlinear evolution of the tearing instability of a cur¬ 
rent sheet (FFE, f = 1, 0 = 1.6) with a = S -1 / 3 . Growth of the 
excited perturbations (the first 8 modes, multiples of ka = 0.03) 
measured on the x average of the B x component. 

4. NONLINEAR SIMULATIONS 

In the present section we investigate the nonlinear 
stages of the evolution of the tearing instability. Since we 
are interested in its late development, where interaction 
and merging of plasmoids is expected, we trigger the in¬ 
stability by selecting an initial spectrum of modes, rather 
than a single one as in the previous set of simulations, 
and we choose a maximum mode number TO max = 10. 
We also choose L y = 1 and S = 10 7 , so the modes with 

/ca ~ 0.029 m; m = 1,10 

are all excited. From the theoretical curves in the pre¬ 
vious section we expect mainly a competition between 
modes to = 3 and m = 4 as the fast growing ones. As 
a reference run, we analyze the instability of a force-free 
equilibrium with constant temperature (FFE, £ = 1), 
and we select the case with /3 = 1.6. This combination 
was shown to provide a linear phase which is the clos¬ 
est to the analytical expectations (see figure [I]). The 
resolution employed for this run is 2048 x 256, which is 
very high if one consider that the code employs high- 
order methods (compact finite-differences along x and 
Fourier transforms along y 1 where periodical boundary 
conditions apply). 

In figure [3] we show the growth in time of the first 8 
excited modes, based on a Fourier analysis along y of the 
B x component averaged in the x direction. Notice that 
towards the end of the linear phase (t ~ 5) the dominant 
modes are m = 3 (ka ~ 0.09), with contributions from 
to = 4, 5, as expected since their growth rates are similar. 
After t ~ 6 — 7 nonlinearities become important starts 
and energy in the longer wavelength modes increases due 
to the merging of smaller wavelength modes: as the m = 
3 mode is the most energetic one at this time, it keeps 
increasing its energy at the expense of higher nr modes. 
At t ~ 9 it has become dominant and the corresponding 
islands start to merge (see the full evolution described 
below): in the spectra this corresponds to a flattening of 
the m = 3 energy curve and to a strong increase first of 


the to = 2 and finally the to = 1 modes (ka = 0.03), the 
largest available island in our y-periodic domain. This 
behavior is compatible with other 2D simulations of the 
tearing instability (see for example figure 3 in Landi et al. 
(2008])). 

The complete evolution is shown in figure |4] where 2D 
snapshots for selected times are provided, of the region 
\x/a\ < 10. Note that the figure has different scales in 
x and y, the x axis is normalized against a, the y axis 
against L , with L/a = S 1//3 ~ 215). Panels on the left 
show the J z current component, whereas panels on the 
right show the distorted field lines superimposed on a 
map of temperature T. It is easy to recognize the end of 
the linear phase of the tearing instability (top row), with 
the dominant m = 3 mode still clearly apparent. 

When the tearing instability growth is over, the nonlin¬ 
ear phase sets in leading to further reconnection events 
and eventually to island coalescence, departing from the 
dominant phase with to = 3. First, due to the attraction 
of current concentrations of the same sign, the X-po ints 


elongate an d stretch along the y direction (Malara et al. 


1991 1992). This process leads to a strong increase of 


the electric current concentrations and thus to the for¬ 
mation of new, elongated current sheets (see the second 
row, details of this phase will be discussed further on). 

Beyond f ~ 9 (third row) the evolution has become 
fully nonlinear and we clearly observe the process leading 
to the creation of a single, large magnetic island as arising 
from coalescence. The situation is very dynamic, espe¬ 
cially near the major X-point where explosive expulsion 
of smaller and smaller islands is observed, typical of the 
plasmoid instability described in Loureiro et al. (2007) 
and subsequent papers (in the spectra, as those shown m 
figure [3j the emergence of plasmoids is revealed by the 
oscillations of modes with to > 1). These islands then 
move towards the largest one, which is continually fed 
and thus further increases its size (approximately with a 
linear behaviour in time, until the boundaries are even¬ 
tually reached) in a sort of inverse cascade eventually 
leading to the largest to = 1 mode. Notice also the tem¬ 
perature enhancement at the reconnection sites. This 
long-term evolution is essentially determined by the im¬ 
posed periodic boundary conditions along y. 

At time t = 9.5 (bottom row) both the current concen¬ 
tration and the plasma temperature around the X-point 
are so high that we need to saturate their values in order 
to retain an appropriate dynamical range in the color 
bar. The initial macroscopic current sheet is basically 
disrupted in a series of highly dynamical features. The 
current and temperature enhancements are stronger at 
the X-point and at the boundaries of the magnetic is¬ 
lands. Moreover, from the major reconnecting site we 
clearly see the production of magnetosonic waves, which 
propagate and soon steepen into shocks. 

At times t > 11 the nonlinearities are so strong that the 
code, which is not conservative, cannot properly resolve 
the shock propagation and the simulation is interrupted. 
In addition, the largest island prevents the formation of 
new, smaller current sheets further feeding it, and the 
subsequent growth appears to proceed slowly, com pati bly 
with the linear growth predicted by [Rutherford (1973). 

Let us now investigate in more detail the"situation right 
after the end of the linear phase of the tearing instability, 
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t= 8.25 



t= 9.50 




-200. 400. 1000. 1600. 0.70 0.77 0.85 0.93 


Fig. 4.— Nonlinear evolution of the tearing instability of a current sheet (FFE, £ = 1, /3 = 1.6) with a = S ~ 1 / : - 5 . On the left panels we 
show the intensity of the J, electric current component, while on the right panels the magnetic fieldlines and the plasma temperature are 
displayed. Notice that the x and y scales are different and that we are zooming the inner region \x\ < 10a, while the computation extends 
out to |x| = 20a. 
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Fig. 5. — Zoom of the most prominent X-point reconnecting region, during the early nonlinear stage of the tearing instability. Colors 
refer to the strength of the J z component of the electric current. On the left panel we show the situation for the reference run with S = 10 7 , 
in the central panel the case S = 10 6 , in the right panel the case S = 10 5 . 


when the islands coalescence is about to set in and the 
local current sheets have just formed and start to further 
evolve. In figure p] (left panel) we show the zoom around 
the X-point (which is just about to develop) near y = 0.2 
at t = 8.25, the first row of figure [4j Here we display the 
electric current by using the same scale in both x and 
y directions, normalized against the macroscopic current 
sheet’s width a = S' -1 / 3 ~ 1/200 ~ 5 x 10 -3 , thus the 
position of the reconnection zone is now expressed as 
y ~ 40a. 

The local current sheet has just formed as the result 
of a stretching process in the y direction (and shrinking 
across the other direction), as a typical output of the 
nonlinear phase of the tearing instability. We are now 
in the phase in which, in the center of the current sheet, 
reconnection is about to take place, leading to a topology 
change in the magnetic structure and to the disruption of 
the current sheet itself, initially into two smaller strips. 
It is very interesting to measure the aspect ratio of this 
reconnection site. From the white box in the figure, de¬ 
fined by the rectangular region where the J z component 
has values which are roughly half those of the central 
peak, we estimate L*/a* — 200, where we have used an 
asterisk to indicate the local values and to differentiate 
with respect to the macroscopic ones (we also recall that 
in the whole paper we identify a as the half width of a 
current sheet). Additional simulations for S = 10 6 and 
S = 10 5 , reported in the second and third panels, lead 
to local aspect ratios of L*/a* ~ 80 and L*/a* ~ 50, 
respectively. 

If we now compare these numbers with S* a , trying to 
find a value of a that fits the data best, it is easy to 
see that the value a = 1/3 is a very good guess. Here 
S* = (L*/L)S is the local Lundquist number, which is 
obviously smaller than the macroscopic one, due to the 
much smaller length of the local current sheet (to be 
measured in each case). Therefore, based on our very 


limited data set, we derive the scaling 
L* /a* = kS* 1/3 , 


( 12 ) 


where k ~ 2.1 — 2.3, that is of order unity, as expected. 

These findings are very important, in our opinion. For 
the first time we clearly see in simulations that, even in 
the nonlinear stages of the tearing instability, the new 
current sheets that form locally become unstable when 
the inverse aspect ratio of these structures reaches the 
critical threshold of a* / L* ~ S'* -1 / 3 , precisely the same 
limit found by PV for the fast reconnection of the initial, 
macroscopic current sheet. After that, a new “ideal” 
tearing instability starts, with time occurring on a faster 
timescale t* a = L* /ca, since typically L* -C L. 

Furthermore, when smaller and smaller scales are pro¬ 
duced nonlinearly, as observed at the time proceeds, 
each time the newly formed local current sheets elongate 
and reach the ir own critical value, that corresponding to 
equation (121, faster and faster reconnection will arise 


producing a cascading, accelerating process: this, we be¬ 
lieve, is the real nature of the plasmoid (or super-tearing) 
instability. 

In order to complete our numerical investigation, we 
have also performed nonlinear simulations for other set¬ 
tings, like the PE initial configuration and /? = 0.1. Some 
changes are seen during the linear stage, because of the 
slightly different growth rates, though the number of 
islands created by the fastest growing mode is always 
m = 3. Coalescence to the m = 2 and m = 1 modes is 
invariably observed, as well as the rapid onset of the plas¬ 
moid instability. The runs with j3 = 0.1 are less robust 
than the reference one with /3 = 1.6, due to the stronger 
compressibility effects and faster growth of shock waves. 
We therefore deem that the nonlinear stages of the fast 
reconnection process described in this section have been 
captured correctly, and the crucial observation that even 
the secondary reconnection events occur when the local 
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current sheets reach an inverse aspect ratio <x as 

for the initial settings, confirms that the physics of the 
“ideal” tearing instability is indeed robust. 

5. CONCLUSIONS 

In the present paper we have studied, by means of 
compressible, resistive MHD simulations, the linear and 
nonlinear stages of the tearing instability of a current 
sheet with inverse aspect ratio a/L ~ S' -1 / 3 , where 
S 1 is the Lundquist number measured on the macro¬ 
scopic scale (the current sheet length) and the asymp¬ 
totic Alfven speed. Our results confirm the linear analy¬ 
sis of PV of the “ideal” tearing mode, while our nonlin¬ 
ear simulations show that the current sheet elongation 
and reconnection follows what appears to be a quasi¬ 
self-similar path, with subsequent collapse, current sheet 
thinning, elongation, destabilization starting from the X- 
points formed in the original sheet. As scales become 
smaller, and the effective Lundquist numbers decrease, 
the dynamical time-scales decrease, leading to explosive 
behavior. 

In particular, we have verified that the secondary 
reconnection events start their fast evolution precisely 
when the critical threshold a/L ~ S' -1 / 3 predicted in 
PV is replicated on the local scale , according to equa¬ 
tion (121. We believe that this could be a universal be¬ 
havior. Though a detailed renormalization-type analysis 
is beyond the scope of this paper, requiring the inclu¬ 
sion of finite inertial-length effects and therefore also the 
presence of guide fields etc., it seems clear that scaling 
phase diagrams for reconnection (as provided for exam¬ 
ple by Cassak & Shay 2012 iCassak & Drake 20131, might 
end up being modified by the presence of this additional 
“ideal” tearing regime and its kinetic generalizations cur¬ 
rently under study. 

Notice that our analysis begins with a current sheet 
which has already reached the critical aspect ratio, but 
it is natural to speculate that because the tearing in¬ 
stability is incredibly slow, at large Lundquist numbers, 
this critical aspect ratio becomes de-facto a trigger for 
explosive energy release. It would be interesting to per¬ 
form large-Reynolds number simulations in which the 
current sheet thins, so that the aspect ratio increases 
with time: in such case the triggering of reconnection 
could be examined in greater detail. Even more realistic 
would be to perform simulatio ns of turbulent reconnec¬ 
tion (see Lazarian et al. 2015 and references therein), 


though it would be hard to reach the accuracy needed to 
recognize our critical scaling in current sheets forming at 
all spatial scales. 

Moreover, here only 2D MHD reconnection has been 
considered, while additional dynamics is expected in 3D 
simulations, with or without the presence of a guide field 
and/or shear flows, where the onset of the secondary 
plasmoid instability may be heavily modified by the in¬ 
teraction with other modes (IQnofri et ah1|2004[ Landi] 


et al. 2008 Bettarini et al. 2009J EancOz^etfarmi 20I2|). 

An important limitation of our simulations is the use 


of on ly a resistive term in Ohm’s law. Tenerani et al. 

2014]) have considered the effects of viscosity on the 
“ideal” tearing scenario, showing that the main result 
is unchanged at Prandtl numbers of order unity, while 
greater critical aspect ratios, even close to the Sweet- 
Parker limit, are possible in more viscous regimes. On 
the other hand, for typical conditions of the solar corona, 
the smaller scales arising in non-linear evolution require 
the inclusion of further kinetic effects, such as the Hall 
term and electron pressure and inertial terms, in Ohm’s 
law. While effects on the linear stability are presently be¬ 
ing considered (Del Sarto et al., submitted) it would be 
important for understanding the partitioning of energy 
between the bulk of the plasma and accelerated particles 
for fully kinetic simulations to be carried out. 

The “ideal” tearing mode occurring on Alfvenic 
timescales naturally applies to the solar coronal plasma, 
where heating is supposed to be caused by the flaring 
activity at all scale s (the Parker nanoflares sce nario) in¬ 
side coronal loops ( Rappazzo et al.|2007 2008). In addi¬ 
tion, the same scenario could most p robably be relevant 


also to relativistic pair plasmas (see Kagan et al. 2015 
and references therein), and in p articular to the model¬ 


ing of Pulsa r Wind Nebulae (e.g. Porth et al.|2014 01mi| 
et al.|[2014) in the light of the recent discovery o f ma.ior 


gamma-ray flares observed in the Crab nebula dTavani 


et al. 


2011[ Komissarov fe Lyutikov 2011 (Jeruttl et al. 


2012[ Buhler fc Blandford|2(jl4). 
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